ST405/ST645 Bayesian Data Analysis

Tutorial Questions (2)

Author

Prof. Niamh Cahill

Solutions

Question 1: Bayesian Model for Multiple Proportions - Email Campaign Click-Through Rates

1. Model Specification

  1. The likelihood for each campaign \(i\) can be modeled as a binomial distribution:

\(\text{y}_i \sim \text{Binomial}(\text{N}_i, \theta_i)\)

where:

  • \(\text{N}_i\) is the number of emails sent for campaign \(i\),

  • \(\text{y}_i\) is the number of emails clicked for campaign \(i\),

  • \(\theta_i\) is the click-through rate (CTR) for campaign \(i\).

  1. We assume a Beta(1,1) prior for each \(\theta_i\), which represents no strong prior belief about the CTR:

\(\theta_i \sim \text{Beta}(1, 1)\)

This is a uniform prior on the interval [0, 1].

2. JAGS Model

  1. Here is the corresponding JAGS code for this model:
CTR_model <- "
model {
  # Loop over the three campaigns
  for (i in 1:n_campaign) {
    # Likelihood: binomial distribution for the number of clicks
    y.i[i] ~ dbin(theta.i[i], N.i[i])
    
    # Prior: Beta(1,1) distribution for the CTR (uniform prior)
    theta.i[i] ~ dbeta(1, 1)
  }
}
"
  1. To implement the JAGS model we need data and paramters to monitor. We can also specify initial values if we wish (JAGS will also do this by default).
  • Data Input for JAGS:
CTR_data <- list(
  y.i = c(15, 20, 12),   # Number of clicks for campaigns A, B, C
  N.i = c(100, 120, 110), # Number of emails sent for campaigns A, B, C
  n_campaign = 3
)
  • Initial Values and Parameters to Monitor:
inits <- function() list(theta.i = runif(3)) # Random initial values for CTRs
params <- c("theta.i") # Monitor the CTRs
  • Putting all of this together and running JAGS
library(rjags)
library(R2jags)

CTR_model <- "
model {
  # Loop over the three campaigns
  for (i in 1:n_campaign) {
    # Likelihood: binomial distribution for the number of clicks
    y.i[i] ~ dbin(theta.i[i], N.i[i])
    
    # Prior: Beta(1,1) distribution for the CTR (uniform prior)
    theta.i[i] ~ dbeta(1, 1)
  }
}
"
CTR_data <- list(
  y.i = c(15, 20, 12),   # Number of clicks for campaigns A, B, C
  N.i = c(100, 120, 110), # Number of emails sent for campaigns A, B, C
  n_campaign = 3
)

inits <- function() list(theta.i = runif(3)) # Random initial values for CTRs
params <- c("theta.i") # Monitor the CTRs

mod <- jags(data = CTR_data, 
            inits = inits, 
            parameters.to.save = params, 
            n.iter = 4000,
            n.burnin = 2000, 
            model.file = textConnection(CTR_model))
Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 3
   Unobserved stochastic nodes: 3
   Total graph size: 11

Initializing model
  1. To check for convergence we can have a quick look at the summary output and trace plots.
mod$BUGSoutput$summary
                 mean         sd        2.5%         25%        50%        75%
deviance   16.2182501 2.39821197 13.49968194 14.48476671 15.6432225 17.2134453
theta.i[1]  0.1574444 0.03536925  0.09512617  0.13208017  0.1556091  0.1798363
theta.i[2]  0.1722803 0.03371329  0.11174441  0.14817670  0.1706998  0.1936532
theta.i[3]  0.1172505 0.02998946  0.06506211  0.09622829  0.1155890  0.1361997
                97.5%     Rhat n.eff
deviance   22.4496301 1.000536  3000
theta.i[1]  0.2338721 1.001354  2300
theta.i[2]  0.2443505 1.000695  3000
theta.i[3]  0.1823926 1.000583  3000
  • Rhat values and ESS look good.
library(coda)
# turn the model into an mcmc object
mod_mcmc <- as.mcmc(mod)
# get trace plot and density
plot(mod_mcmc)

3. Interpretation of Results

  • After running the MCMC sampling (e.g., using 4,000 iterations and discarding the first 2,000 as burn-in), you will obtain samples from the posterior distributions of \(\theta_i\).

  • (a) Plot the posterior distributions of the CTRs for the three campaigns using a package like tidybayes in R (or coda or the base plotting functions).

library(tidybayes)
library(tidyverse)
## get the output in matrix format
mcmc.matrix <- mod$BUGSoutput$sims.matrix
## get indexes for theta
campaign_index <- 1:3

## plots
mcmc.matrix %>%
  spread_draws(theta.i[campaign_index]) %>% 
  ggplot(aes(y = factor(campaign_index, labels = c("A","B", "C")), x = theta.i)) +
  stat_halfeye() +
  ylab("Campaign")

  • (b) Compare the distributions to assess the relative effectiveness of the campaigns. It appears that the posterior distribution for \(\theta_B\) is shifted to the right compared to \(\theta_A\) and \(\theta_C\), suggesting that Campaign B is likely more effective, i.e., the click through rate is estimated to be higher. Note however, there is an overlap in terms of uncertainty in the estimated click through rate.

4. Additional Analysis

  • (a) Calculate the posterior probability that Campaign B has a higher CTR than Campaign A and Campaign C:

\(P(\theta_B > \theta_A) = \frac{1}{S} \sum_{s=1}^{S} I(\theta_B^{(s)} > \theta_A^{(s)})\)

where \(S\) is the number of posterior samples, and \(I(\cdot)\) is an indicator function. You can do similar calculations for \(P(\theta_B > \theta_C)\).

We can implement this in R as follows

theta_samps <- mcmc.matrix %>%
                spread_draws(theta.i[campaign_index]) 

theta_A <- theta_samps %>% filter(campaign_index == 1) %>% pull(theta.i)
theta_B <- theta_samps %>% filter(campaign_index == 2) %>% pull(theta.i)
theta_C <- theta_samps %>% filter(campaign_index == 3) %>% pull(theta.i)

sum(theta_B > theta_A)/length(theta_B)
[1] 0.6276667
sum(theta_B > theta_C)/length(theta_B)
[1] 0.8923333

Interpretation:

  • The posterior probability that \(\theta_B > \theta_A\) is 0.61, indicating that there is a 61% chance Campaign B has a higher CTR than Campaign A.
  • The posterior probability that \(\theta_B > \theta_C\) is 0.89, indicating that there is a 89% chance Campaign B has a higher CTR than Campaign C.